'''
注意本方法中的数据为直接的遥感影像的数据
这里尝试结合图分割与SLIC的变形结合
前提假设:
1.建筑物的边界应是闭合的。--> 边缘本身应该是的闭合，但是在实际中发现，本应是闭合的边缘实际中的，并不存在的闭合，而是开放性的，因此这一条假设释放掉。
2.神经网络的FCN计算结果可以等价于像素隶属于某一类的模糊度量值（值范围：0~1）
3.任何相邻两个邻域部件之间被边界隔绝，即边界是不可跨越的 
4.边界处的曲率不同于部件内部的曲率值
5.
图表达结合传统的矩阵形式
step 1: 使用的基于曲率的变形canny算法，提取边缘
step 2: 基于最小树生长算法，构建子图簇。度量各个子图的特性，并计算隶属度。

使用传统的边缘提取算子，强化图像中边缘性，并且利用建筑物边界的闭合性的滤除一部分非边缘边界。
step 1. 边缘边界检测：
    1. canny 常规算子
    2. canny 曲率算子
step 2. 非边缘边界检测：
    1. 其中边界的闭合性，使用图论中的闭合环检测。
        a. 直接删除 度为1 的节点
step 3. 最小生成树：
    1. 图的结构定义
        a. 像素表示单个顶点
        b. 边表示像素与子团之间的相似性
    2. 合并规则：
        a. 边界闭合性作为的空间约束
        b. 像素的特征计算
        c. 深度学习结果作为隶属度函数的度量
        d. 像素相似性计算
        e. 像素与临近团的相似性计算
    3. 分裂规则：
        a. 度量每个团簇的熵值
        b. 分裂方法
    4. 最小生成树终止标准
        a. 以整体熵值最小作为的度量标准
        b. 所有的像素全部已经调整到了合适的子团中
step 4. 图像中的语义隶属度分割
    1. 度量图元的语义隶属度
        a. 平均值
        b. 值域范围（最小值，最大值）
        c. 标准差
        d. 整体概率分布（从大到小的序列排序）
    2. 最近邻合并同一类型子团
        a. 遵循联通结构
'''

import numpy as np
from numpy.core.defchararray import array
from numpy.ma.core import set_fill_value
from scipy import sparse
from scipy import stats # 用于拟合gamma分布
from scipy.sparse import csr_matrix,lil_matrix,csc_matrix,coo_matrix
from scipy.sparse.construct import random # 稀疏矩阵
import torch
from torch._C import Graph
import torch.nn.functional as F
import math
import scipy.io as scio
from matplotlib import pyplot as plt
import os
import random
from osgeo import gdal
import copy

class NetPredict(object):
    '''
    这个方法用来加载与处理网络模型的预测结果
    '''
    def __init__(self):
        super().__init__()
        pass
    def LoadModel(self,modelName,modelcpk_path,num_cls,input_channels):
            '''
            加载模型
            args:
                modelName 需要参与计算的模型的名称
                modelcpk_path 训练好模型的保存路径
            return:
                model : 加载好对应权重的模型
            '''
            model=get_models(modelName,num_class=num_cls,channels_num=input_channels)
            epoch,model,best_iou=load_ckpt(model,None,modelcpk_path,'cuda')
            self.model=model
            self.epoch=epoch
            self.best_iou=best_iou
            return 1
    def predict(self,img):
        '''
        利用预加载好的模型，预测图像
        '''
        return self.model(img)




class GraphBase(object):
    '''
    定义基本的图网络操作方法，与基本结构
    '''
    
    def __init__(self,node_data,node_id_arr,node_id_xy,edges):
        super().__init__()
        self.node_data=node_data
        self.node_id_arr=node_id_arr
        self.node_id_xy=node_id_xy
        self.edges=edges
        self.node_num=node_id_arr.shape[0]*node_id_arr.shape[1] # 统计点数

    def GetNodeData(self,nid):
        xy=self.getNode_XY(nid)

    def FindNode_Near(self,node_id):
        '''
        查找到当前节点的周围邻接节点 Start -> End
        '''
        return  self.edges[node_id,:].nonzero()[1].tolist()# 获得所有相连的边

    def FindNode_End(self,node_id):
        '''
        查找尾在当前节点的 节点 End-> Start
        '''
        return self.edges[:,node_id].nonzero()[0].tolist() 

    def getNodes(self):
        '''
        将当前节点序列化
        '''
        return np.copy(self.node_id_arr).reshape(-1).tolist()
    def getNode_XY(self,nid):
        '''
        获得当前坐标的xy
        '''
        return self.node_id_xy[nid,:].tolist()

    def GenerateNetwork(img_data):
        '''
        创建对应的图形结构
        使用稀疏矩阵构建对应的结构，
        网络图结构定义为dict，而不是矩阵结构，
        G=(N,E) N 表示节点，E 表示边
        其中N用矩阵 表达+ 矩阵坐标

        E用dict表达+稀疏矩阵
        '''
        h,w=img_data.shape
        node_id_arr=np.array(list(range(h*w))).reshape(h,w).astype(np.int) # Node id 的节点板
        w_arr=np.ones((h,w),dtype=np.int)*np.array(list(range(w))) # x
        h_arr=np.ones((w,h),dtype=np.int)*np.array(list(range(h))) # y
        h_arr=h_arr.transpose()
        w_arr=w_arr.reshape(-1)
        h_arr=h_arr.reshape(-1)
        h_w_arr=np.array([h_arr,w_arr]).transpose().astype(np.int)
        # 记录id之间的连接关系,每一个像素有8 邻域
        edge_dict={}
        node_id_pad_arr=np.pad(node_id_arr,pad_width=1,mode='edge') # 边界外扩展
        # 构建8邻域的边界点
        LU_node_list=node_id_pad_arr[:-2,:-2].reshape(-1).astype(np.int).tolist() #  -1,-1
        MU_node_list=node_id_pad_arr[:-2,1:-1].reshape(-1).astype(np.int).tolist()# -1,0
        RU_node_list=node_id_pad_arr[:-2,2:].reshape(-1).astype(np.int).tolist()#   -1,1

        LM_node_list=node_id_pad_arr[1:-1,:-2].reshape(-1).astype(np.int).tolist() # 0-1
        MM_node_list=node_id_pad_arr[1:-1,1:-1].reshape(-1).astype(np.int).tolist()# 0,0
        RM_node_list=node_id_pad_arr[1:-1,2:].reshape(-1).astype(np.int).tolist()  # 0,1

        LD_node_list=node_id_pad_arr[2:,:-2].reshape(-1).astype(np.int).tolist() #   1,-1
        MD_node_list=node_id_pad_arr[2:,1:-1].reshape(-1).astype(np.int).tolist() #  1,0 
        RD_node_list=node_id_pad_arr[2:,2:].reshape(-1).astype(np.int).tolist() #    1,1

        # 初始边界值
        weight_a_b_list=np.ones((h,w)).reshape(-1).tolist()

        # 构建稀疏阵
        node_len=len(MM_node_list)
        edge_coo=coo_matrix((weight_a_b_list,(MM_node_list,LU_node_list)),shape=(node_len,node_len))
        edge_coo=coo_matrix((weight_a_b_list,(MM_node_list,MU_node_list)),shape=(node_len,node_len))+edge_coo
        edge_coo=coo_matrix((weight_a_b_list,(MM_node_list,RU_node_list)),shape=(node_len,node_len))+edge_coo
        edge_coo=coo_matrix((weight_a_b_list,(MM_node_list,LM_node_list)),shape=(node_len,node_len))+edge_coo
        edge_coo=coo_matrix((weight_a_b_list,(MM_node_list,MU_node_list)),shape=(node_len,node_len))+edge_coo
        edge_coo=coo_matrix((weight_a_b_list,(MM_node_list,RM_node_list)),shape=(node_len,node_len))+edge_coo
        edge_coo=coo_matrix((weight_a_b_list,(MM_node_list,LD_node_list)),shape=(node_len,node_len))+edge_coo
        edge_coo=coo_matrix((weight_a_b_list,(MM_node_list,MD_node_list)),shape=(node_len,node_len))+edge_coo
        edge_coo=coo_matrix((weight_a_b_list,(MM_node_list,RD_node_list)),shape=(node_len,node_len))+edge_coo
        edge_coo=edge_coo+edge_coo.T
        edge_coo=(edge_coo>0)*1
        edge_csr=edge_coo.tocsr()
        # edge_coo 结果为对称阵。
        return  GraphBase(img_data,node_id_arr,h_w_arr,edge_csr)


class GraphCUT(object):
    '''
    图分割算法基本理论架构基本方法结构。
    注意这里因为需要考虑 子图 计算的因素。
    特意说明图的结构如下:
        1. 节点： Node :id_xy,id_arr
        2. 边界： start_edge ：{start:[end]} A->B
                 end_edge:{end:[start]} B->A
        

    算法粗略结构图
    step 1 创建图
    step 2 计算节点与周围节点的相近性
    step 3 合并最大可能的节点簇
    step 4 计算整体簇的熵值
    step 5 当熵值最小时结束
    '''
    def __init__(self):
        '''
        图割法初始化

        args:
            img_data :np.darray h,w 图像，用于生成点集。
        '''
        super(GraphCUT,self).__init__()
        pass
    

    def computerEntroy(self,node_ls):
        '''
        计算当前节点列表中的整体熵值情况
        '''
        return 0
        pass

    def InitSubGraph(self):
        '''
        初始化子图，这里主要确定节点与子团之间的关联
        '''
        self.sub_graph={} # 定义子图
        for i in self.graph.getNodes():
            [x,y]=self.graph.getNode_XY(i)
            self.sub_graph[int(i)]={'nodes':[int(i)],'entropy':-1,'predict':self.NodeData['img_predict'][x,y]} # 其中的 entropy 表示当前子图的熵，predict 为隶属度差值平均值 
            # 度量熵值
            self.sub_graph[int(i)]['entropy']=self.computerEntroy(self.sub_graph[int(i)]['nodes'])

        # 定义子图膜版
        self.sub_graph_arr=np.copy(self.graph.node_id_arr) # 初始化子图归属 根据对应的节点位置，直接确定子图情况
        self.node_subgraph={} # 判断单个节点邻接性

    def GetSubGraphID(self,nid):
        '''
        返回当前节点id 对应的子图
        '''
        n_h,n_w=self.graph.getNode_XY(nid)
        return self.sub_graph_arr[n_h,n_w] # 返回对应子图的id

    def GetNearSubGraph(self,nid):
        '''
        获取当前节点周围的所有的子图,并计算相似性
        '''
        [N_h,N_w]=self.graph.getNode_XY(nid) # 
        Near_nids=self.graph.FindNode_Near(nid) # 获取当前邻接点
        # 求解对应的子图
        result=[]
        for nid_near in Near_nids:
            [n_h,n_w]=self.graph.getNode_XY(nid_near)
            result.append(self.sub_graph_arr[n_h,n_w]) # 追加子图
        result=list(set(result)) # 
        return result

    def ComputerSimilarityNode(self,A_id,B_id,megerEdge=True):
        '''
        计算两个节点之间的相似性，
        '''
        # 度量两个节点之间 光谱的相似性---这里先使用相关性系数，
        [A_h,A_w]=self.graph.getNode_XY(A_id)
        [B_h,B_w]=self.graph.getNode_XY(B_id)
        # 判断 B_id 是否是边界，若是边界，直接的返回空置，表示当前两个节点不可连接
        if self.NodeData['img_edge'][B_h,B_w]==1 and megerEdge: # 节点不参与计算
            return None
        
        A_ori_data=self.NodeData['img_data'][A_h,A_w,:]
        B_ori_data=self.NodeData['img_data'][B_h,B_w,:]
        [ori_distance,ori_confidence]=np.corrcoef(A_ori_data,B_ori_data)[0,:].tolist() # 相关系数返回  相关系 与 置信度
        # 这里只考虑正相关的情况 
        if ori_confidence<0.9 or ori_confidence <0.7:  # --------------------------------------------------------------------------------------这里存在0.9 与 0.7 为超参数
            ori_confidence=0 # 相关系数为0
        # 度量各个波段性的相似性
        A_feature=self.NodeData['img_process'][A_h,A_w,:]
        B_feature=self.NodeData['img_process'][B_h,B_w,:]
        [feature_distance,feature_confidence]=np.corrcoef(A_feature,B_feature)[0,:].tolist()
        if feature_confidence<0.9 and feature_distance<0.7: #---------------------------------------------------------------------------------这里0.9 与 0.7 同样为超参数
            feature_distance=0
        
        # 度量隶属度
        A_predict=self.NodeData['img_predict'][A_h,A_w]
        B_predict=self.NodeData['img_predict'][B_h,B_w]
        predict_distance=abs(A_predict-B_predict) # 度量两者差别
        predict_distance=1-predict_distance # 归一化 保证趋近1 与前两个同样
        # 同一范围约束
        ori_distance=1-ori_distance
        feature_distance=1-feature_distance
        predict_distance=1-predict_distance

        return ori_distance,feature_distance,predict_distance # 所有的度量范围 [0,1] 数量越小，越相近

    def ComputerSimilarityNodeAndSubgraph(self,nid,subgraph_id):
        '''
        计算当前节点与相邻子团的相似性
        '''
        ori_distance=[] # 原始波段相关性，越小越好
        predict_distance=[] # 预测波段差异性 越小越好
        feature_distance=[] # 预测特征相关性 越小越好
        predict_mean=-1 # 预测平均值与节点预测值的差别
        Entropy_delta=0
        for Bid in self.sub_graph[subgraph_id]['nodes']:
            # 整体度量值
            if Bid==nid: # 计算除了本体之外的结果
                continue 
            ori_distance_temp,feature_distance_temp,predict_distance_temp=self.ComputerSimilarityNode(nid,Bid)
            ori_distance.append(ori_distance_temp)
            feature_distance.append(feature_distance_temp)
            predict_distance.append(predict_distance)
            pass
        predict_mean=abs(self.NodeData['img_predict']-self.sub_graph[subgraph_id]['predict']) # 平均值之差
        # 计算考虑nid 融入到subgraph 之中，熵值的变化。

        if nid in self.sub_graph[subgraph_id]['nodes']: # 此时不需要变更熵值
            Entropy_delta=0
        else:
            # 变更熵值由两部分组成【离开子团时，熵值减少。增加子团时，熵值增加】
            # 获得当前子节点
            n_sub_graph_id=self.GetSubGraphID(nid)
            nodes_d=[node_id_in_same_subgraph   for node_id_in_same_subgraph in self.sub_graph[n_sub_graph_id]['nodes']  if node_id_in_same_subgraph!=nid ]
            nodes_add=[node_id_in_near_subgraph for node_id_in_near_subgraph in self.sub_graph[subgraph_id]['nodes']]
            nodes_add.append(nid)
            entropy_add=self.computerEntroy(nodes_add)
            # 判断当前网络的结构情况
            if entropy_add>1 : #--------这里的1是一个超参数，用于控制当前超像素团的聚集情况-------------------------------------------------------------------------------------------------------------------------------------------------------------- 
                return None

            # 计算熵值变化
            Entropy_delta=self.sub_graph[n_sub_graph_id]['entropy']-self.computerEntroy(nodes_d)+self.sub_graph[subgraph_id]['entropy']-entropy_add # 计算熵值变化
        # 考量三个方向之差 与 熵值变化
        result=0 # 最终的结果------------------------------------------------------------------------------------------------------------------------------------------------------
        
        return result

    def ComputerSimilarity(self):
        '''
        计算像素与相邻子团的相似性,并返回最大的相邻子团以及对应的距离值
        '''
        # 记录列表,用于构建 相似性矩阵
        nodes_subgraph_data={}
        graph_nodes_list=self.graph.getNodes()
        for nid in graph_nodes_list:
            [n_h,n_w]=self.graph.getNode_XY(nid) # 确定当前节点的坐标
            sub_graph_ls=self.GetNearSubGraph(nid) # sub_graph_ls 获得所有的临近相似的度
            max_distanct=math.inf
            max_subgraph_id=0
            for sub_graph_id in sub_graph_ls:
                temp=self.ComputerSimilarityNodeAndSubgraph(nid,sub_graph_id)
                if temp is None:
                    continue
                if temp<max_distanct:
                    max_subgraph_id=sub_graph_id
            nodes_subgraph_data[(nid,max_subgraph_id)]=copy.deepcopy(temp)
        return nodes_subgraph_data

    def UpdataSubgraph(self,nodes_subgraph_data):
        for k in nodes_subgraph_data:
            # 重新平移
            (nid,subgraph_id)=k
            n_subgraph_id=self.GetSubGraphID(nid)
            self.sub_graph[n_subgraph_id]['nodes'].remove(int(nid))
            self.sub_graph[subgraph_id]['node'].append(nid)
            n_h,n_w=self.graph.getNodes(nid)
            self.sub_graph_arr[n_h,n_w]=subgraph_id
        # 分别更新节点结果
        for k in self.sub_graph:
            if len(self.sub_graph[k]['nodes'])==0:
                self.sub_graph[k]['predict']=0
                self.sub_graph[k]['entropy']=0
                continue
            mean_predict=0
            for nid in self.sub_graph[k]['nodes']:
                n_h,n_w=self.graph.getNode_XY(nid)
                mean_predict=mean_predict+self.NodeData['img_predict'][n_h,n_w]
            mean_predict=mean_predict/len(self.sub_graph[k]['nodes'])
            self.sub_graph[k]['predict']=mean_predict
            self.sub_graph[k]['entropy']=self.computerEntroy(self.sub_graph[k]['nodes']) # 计算当前节点的熵
        
    def ComputerMeanEntropy(self):
        '''
        计算整体熵值的平均值
        '''
        sum_entropy=0
        d=0
        for k in self.sub_graph:
            sum_entropy=sum_entropy+self.sub_graph[k]['entropy']
            d=d+1
        mean_entropy=sum_entropy/d
        return mean_entropy

    def computer(self,img_data,h,w):
        '''
        图分割计算方法
        img_data:dict  图像数据 数据结构{img_data,img_predict,img_edge,img_preprocess}
        h,w：图像大小
        '''
        # 构建网络
        self.graph=GraphBase.GenerateNetwork(np.zeros((h,w),dtype=np.int)) # 创建图临近关系
        self.NodeData=img_data # 定义数据结构
        self.InitSubGraph() # 初始化子图
        '''图分割法开始构建'''
        last_entropy=10000
        while epoch<100:
            # step 1 计算点-子图关系结构
            nodes_subgraph_data=self.ComputerSimilarity() # 计算节点与各个子图的相似度
            self.UpdataSubgraph(nodes_subgraph_data) # 更新各个子图
            mean_entroyp=self.ComputerMeanEntropy() # 计算整体的平均熵值
            if abs(mean_entroyp-last_entropy)<1:
                break
            last_entropy=mean_entroyp # 计算整体熵值变化
            epoch=epoch+1
            pass
        
class canny(object):
    '''
    算法基本步骤：
    step 1: 图像灰度图化（但是本文所涉及到的影像数据本身就是灰度值，跳过了这一步）
    step 2: 图像滤波平滑图像（这里使用了高斯曲率滤波，高斯平滑滤波，建议高斯曲率平滑滤波算子）
    step 3: 图像梯度值与方向计算，获取梯度值与梯度方向（这里使用sobel微分差分算子 作为梯度检测与梯度方向计算）
    step 4: 图像非极大值抑制
    step 5: 阈值滤波这里使用了的OTSU法查找 二分类阈值
    step 6: 滞后边界跟踪 ---这里采用最小生成树算法，
    step 7: 返回结果
    '''
    def __init__(self,pred_process,smooth_flter="wmc",edge_filter="sobel"):
        super(canny,self).__init__()
        self.pred_process=pred_process # 图像预处理

        self.smooth_filter=canny.wmcfilter  if smooth_flter=="wmc" else canny.gaussfilter # 曲率平滑
        self.edge_filter=canny.sobelfilter if edge_filter=="sobel" else canny.linearCurvaturefilter # 边界求解
        self.NMS=canny.NMS

    # 平滑滤波 ---------------------------------
    def wmcfilter(img_gray):
        '''
        曲率滤波的线性化运算
        { 
        @article{GONG2019,
        title = "Weighted mean curvature",
        journal = "Signal Processing",
        volume = "164",
        pages = "329 - 339",
        year = "2019",
        issn = "0165-1684",
        doi = "https://doi.org/10.1016/j.sigpro.2019.06.020",
        url = "http://www.sciencedirect.com/science/article/pii/S0165168419302282",
        author = "Yuanhao Gong and Orcun Goksel",}
        }
        '''
        ori_img_gray=np.copy(img_gray)
        # 定义卷积核
        h,w=img_gray.shape
        h1=np.array([[1,1,0],[2,-6,0],[1,1,0]])/6.0
        h2=np.copy(np.transpose(h1))
        h3=np.copy(np.fliplr(h1))
        h4=np.copy(np.flipud(h2))
        h5=np.array([[2,4,2],[4,-12,0],[1,0,0]])/12.0
        h6=np.copy(np.fliplr(h5))
        h7=np.copy(np.flipud(h5))
        h8=np.copy(np.fliplr(h7))
        # 构建kernel
        kernel=torch.from_numpy(np.array([h1,h2,h3,h4,h5,h6,h7,h8]).reshape(8,1,3,3)).float()
        dist=np.zeros((8,1,h,w))
        # 使用边缘填充
        dist=torch.from_numpy(dist)
        img_gray=torch.from_numpy(img_gray).reshape(1,1,h,w).float()
        img_gray=torch.nn.ReplicationPad2d(1)(img_gray)
        dist=F.conv2d(img_gray,kernel,padding=0) 
        dist=dist.numpy().reshape(8,h,w) # 1,8,h,w
        dist=dist.reshape(8,h*w) # 修正为 8,h*w 方便计算列向量
        min_dist_abs=np.min(np.abs(dist),axis=0) # 获得距离的绝对值
        abs_mask=dist=min_dist_abs
        abs_mask=np.sum(abs_mask,axis=0)==0 # 表示值为负数
        abs_mask=-2*abs_mask
        abs_mask=abs_mask+1
        min_dist_abs=min_dist_abs*abs_mask
        min_dist_abs=min_dist_abs.reshape(h,w)
        ori_img_gray=ori_img_gray+min_dist_abs
        return ori_img_gray

    def gaussfilter(img_gray,sigmoid=1.5,kernel_size=3):
        '''
        高斯滤波
        args:
        img_gray    :   灰度图像
        sigmoid :   标准差值
        kernel_size :   核的大小
        '''
        h,w=img_gray.shape
        # 生成mask，kernel --对称核
        kernel=np.zeros((kernel_size,kernel_size))
        edge_size=int((kernel_size-1)/2)
        for i in range(0,int(edge_size+1)):
            for j in range(0,int(edge_size+1)):
                kernel[int(i+edge_size),int(j+edge_size)]=math.exp(-1*(i**2+j**2)/(2*(sigmoid**2)))/(2*math.pi*(sigmoid**2)) #计算对应的高斯值      
        kernel[:edge_size,:]=np.flipud(kernel)[:edge_size,:]
        kernel[:,:edge_size]=np.fliplr(kernel)[:,:edge_size]
        kernel=kernel/np.sum(kernel) # 归一化  保证各权重之和为 1
        kernel=torch.from_numpy(kernel).reshape(1,1,kernel_size,kernel_size).float()
        img_gray=torch.from_numpy(img_gray).reshape(1,1,h,w).float()
        # 
        img_gray=torch.nn.ReplicationPad2d(1)(img_gray)
        filter_result=F.conv2d(img_gray,kernel,padding=0)
        filter_result=filter_result.reshape(h,w)
        filter_result=filter_result.numpy()
        return filter_result

    # 边缘滤波-------------------------
    def linearCurvaturefilter(img_gray):
        '''
        此方法用于求解平均曲率，来源与 wmcfilter 论文作者的博士论文中的6.12 公式
        '''
        kernel=np.array([[-1,5,-1],[5,-16,5],[-1,5,-1]])/16.0
        kernel=torch.from_numpy(kernel).float().reshape(1,1,3,3) # 卷积核
        h,w=img_gray.shape
        img_gray=torch.from_numpy(img_gray).reshape(1,1,h,w).float()
        # padding
        img_gray=torch.nn.ReplicationPad2d(1)(img_gray)
        img_gray=F.conv2d(img_gray,kernel,padding=0) 
        img_gray=img_gray.numpy().reshape(h,w)
        return img_gray

    def sobelfilter(img_gray):
        '''定义sobel滤波初始算子'''
        kernel_x=np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
        kernel_y=np.array([[1,2,1],[0,0,0],[-1,-2,-1]])
        kernel_x=torch.from_numpy(kernel_x).reshape(1,1,3,3).float()
        kernel_y=torch.from_numpy(kernel_y).reshape(1,1,3,3).float()
        # 考虑数据运算
        h,w=img_gray.shape
        img_gray=torch.from_numpy(img_gray).reshape(1,1,h,w).float()
        img_gray_ex=torch.nn.ReplicationPad2d(1)(img_gray)
        G_x=F.conv2d(img_gray_ex,kernel_x,padding=0)
        G_y=F.conv2d(img_gray_ex,kernel_y,padding=0)
        G=torch.pow(G_x,2)+torch.pow(G_y,2)
        G=torch.pow(G,0.5) # 求解梯度
        print(torch.max(G),torch.min(G))
        # 为避免div 0 这里进行数据规划
        G_mask=G>0 # 选择所有大于0 的值
        G=(G==0)+G # 所有大于0 的值+1 
        
        G_theta=G_y/G # 求解sin值，且最小值为0 
        G_theta=G_theta*G_mask # 将所有梯度值为0的值，置为0

        print(torch.max(G_theta),torch.min(G_theta))
        G_theta=torch.asin(G_theta) # 求解梯度方向 ,弧度制表达
        
        G=G.numpy().reshape(h,w)
        G_theta=G_theta.numpy().reshape(h,w)
        return [G,G_theta] # 返回梯度值与梯度方向
    def NMS(G,G_theta):
        '''
        非极大值抑制，找寻像素点局部最大值，将非极大值点所对应的灰度值置0
        G： 梯度值
        G_theta: 梯度方向（弧度制）
        算法基本流程：
        Step 1 选中的x,y
        Step 2 选中的周围采样点 g1,g2,gn
        step 3 求解相交点 ptemp1 ptemp2
        step 4 求解相交点 pt1,pt2 的双线性差值
        step 5 比较 g_xy pt1 pt2 的值，如果g_xy 是最大值，则保留，不是则置为0
        step 6 重复过程 ，直至所有的可能性检索完成
        注：
        这里使用循环，因此可能整体性能会下降的巨大，因此这里改用基于torch的矩阵计算型，避免循环的使用
        '''
        h,w=G.shape
        G=torch.from_numpy(G).float() # 梯度值 转换成 tensor
        G_theta=torch.from_numpy(G_theta).float() # 梯度方向值 转换为tensor
        # 填充边界
        G_Ex=torch.nn.ReplicationPad2d(1)(G.reshape(1,1,h,w)).reshape(h+2,w+2) # 扩充梯度值边界
        G_theta_Ex=torch.nn.ReplicationPad2d(1)(G_theta.reshape(1,1,h,w)).reshape(h+2,w+2) # 扩充梯度方向值边界
        # 构建偏移矩阵
        offset_x=np.ones((h,w))*np.array(list(range(1,int(w+1)))) # 构建偏移矩阵 x 
        offset_y=np.ones((w,h))*np.array(list(range(1,int(h+1)))).transpose() # 构建偏移矩阵 y 
        offset_x=torch.from_numpy(offset_x).float() # 转换为 tensor
        offset_y=torch.from_numpy(offset_y).float() # 转换为 tensor
        # 求解梯度方向值，并构建梯度方向差值坐标 pt1,pt2。 其中pt2与pt1中心对称，假定中心为（0，0）,求解对应的插值情况
        # 分成四种情况（-pi/2 ,-pi/4],(-pi/4,pi/4],(pi/4,pi/2)，{-pi/2,pi/2} ,
        # 其中pi/2与-pi/2 在几何意义上一致，只讨论pi/2 的情况
        pt1_x=G_theta.clone()*0
        pt1_y=G_theta.clone()*0
        tan_theta=torch.tan(G_theta) # 计算所有角度的tan值，方便作为后续进行计算
        tan_theta_0=(tan_theta==0)*1+tan_theta
        tan_theta_reciprocal=1/tan_theta_0 # 获得所有零的倒数
        # 计算（-pi/2,-pi/4） 此时y=-1 ，x=-1/tan(theta)
        mask_theta=1.0*(G_theta<(math.pi/-4))*(G_theta>(math.pi/-2)) # 条件约束版
        temp_x= -1*mask_theta*tan_theta_reciprocal # 这里表示第四象限
        pt1_x=pt1_x+temp_x # -1/tan(theta)
        pt1_y=pt1_y+mask_theta*-1 # -1 
        # 计算 （-pi/4,pi/4） 此时x=1,y=tan(theta) 
        mask_theta=1.0*(torch.abs(G_theta)<math.pi/4) # mask 约束
        temp_y=mask_theta*tan_theta
        pt1_x=pt1_x+mask_theta # 1
        pt1_y=pt1_y+temp_y # tan(theta)
        # 计算 （pi/4,pi/2） 此时y=1,x=-1/tan(theta)
        mask_theta=1.0*(G_theta>math.pi/4)
        temp_x=mask_theta*tan_theta_reciprocal
        pt1_x=pt1_x+temp_x # 1/tan(theta)
        pt1_y=pt1_y+mask_theta # 1
        # 单独计算G_theta=pi/2 or -pi/2 的情况
        mask_theta=(torch.abs(G_theta)==math.pi/2)*1.0 # 获得所有的特殊值
        mask_theta_no_pi_div_2=(torch.abs(G_theta)<math.pi/2)*1.0 # 其余值
        pt1_x=pt1_x*mask_theta_no_pi_div_2 # 置0
        pt1_y=pt1_y*mask_theta_no_pi_div_2 # 置0
        pt1_x=pt1_x # 此时 x=0
        pt1_y=pt1_y+mask_theta # 此时 y=1
        # 求解对称阵 pt2
        pt2_x=pt1_x*-1 # x的对称
        pt2_y=pt1_y*-1 # y的对称
        # 偏差修正 [dx,dy]+[x,y]=[x+dx,y+dy] 偏差方程
        pt1_x=pt1_x+offset_x 
        pt1_y=pt1_y+offset_y
        pt2_x=pt2_x+offset_x
        pt2_y=pt2_y+offset_y
        # 偏差点 pt1
        # 偏差点 pt2
        # 求解采样点 pt1 的插值结果
        # 针对 pt1 这里采取使用的双线性插值
        # 双线性插值公式 
        # f(p)= b*d*Q11+a*d*Q21+b*c*Q12+a*c*Q22
        # a=x-x1
        # b=x2-x
        # c=y-y1
        # d=y2-y
        # Q11=(x1,y1)
        # Q12=(x1,y2)
        # Q21=(x2,y1)
        # Q22=(x2,y2)
        pt_x1=torch.floor(pt1_x) # 向下取整 获得 x1
        pt_x2=torch.ceil(pt1_x)  # 向上取整 获得 x2
        pt_y1=torch.floor(pt1_y) # 向下取整 获得 y1
        pt_y2=torch.ceil(pt1_y)  # 向上取整 获得 y2 
        # 构建计算系数 a，b,c,d
        a=pt1_x-pt_x1 # (x-x1)
        b=pt_x2-pt1_x # (x2-x)
        c=pt1_y-pt_y1 # (y-y1)
        d=pt_y2-pt1_y # (y2-y)
        a=a.reshape(-1) 
        b=b.reshape(-1)
        c=c.reshape(-1)
        d=d.reshape(-1)
        # 构建采样坐标
        pt_x1=pt_x1.numpy().astype(np.int32).reshape(-1) # x1 
        pt_x2=pt_x2.numpy().astype(np.int32).reshape(-1) # x2
        pt_y1=pt_y1.numpy().astype(np.int32).reshape(-1) # y1
        pt_y2=pt_y2.numpy().astype(np.int32).reshape(-1) # y2
        # 整体采样
        Q11=G_Ex[pt_x1,pt_y1].reshape(-1)
        Q12=G_Ex[pt_x1,pt_y2].reshape(-1)
        Q21=G_Ex[pt_x2,pt_y1].reshape(-1)
        Q22=G_Ex[pt_x2,pt_y2].reshape(-1)
        # 开始计算
        Fp1=b*d*Q11+a*d*Q21+b*c*Q12+a*c*Q22
        Fp1=Fp1.reshape(h,w)
        # 求解采样点 pt2 的插值结果 双线性插值--------------------------
        pt_x1=torch.floor(pt2_x) # 向下取整 获得 x1
        pt_x2=torch.ceil(pt2_x)  # 向上取整 获得 x2
        pt_y1=torch.floor(pt2_y) # 向下取整 获得 y1
        pt_y2=torch.ceil(pt2_y)  # 向上取整 获得 y2 
        # 构建计算系数 a，b,c,d
        a=pt2_x-pt_x1 # (x-x1)
        b=pt_x2-pt2_x # (x2-x)
        c=pt2_y-pt_y1 # (y-y1)
        d=pt_y2-pt2_y# (y2-y)
        a=a.reshape(-1) 
        b=b.reshape(-1)
        c=c.reshape(-1)
        d=d.reshape(-1)
        # 构建采样坐标
        pt_x1=pt_x1.numpy().astype(np.int32).reshape(-1) # x1 
        pt_x2=pt_x2.numpy().astype(np.int32).reshape(-1) # x2
        pt_y1=pt_y1.numpy().astype(np.int32).reshape(-1) # y1
        pt_y2=pt_y2.numpy().astype(np.int32).reshape(-1) # y2
        # 整体采样
        Q11=G_Ex[pt_x1,pt_y1].reshape(-1)
        Q12=G_Ex[pt_x1,pt_y2].reshape(-1)
        Q21=G_Ex[pt_x2,pt_y1].reshape(-1)
        Q22=G_Ex[pt_x2,pt_y2].reshape(-1)
        # 开始计算
        Fp2=b*d*Q11+a*d*Q21+b*c*Q12+a*c*Q22
        Fp2=Fp2.reshape(h,w)    
        # Fp1,Fp2 插值点，用于非极大值比较
        # 当p>Fp1 && p> Fp2 保留，其余置零
        Fp1=(Fp1<G)*1.0
        Fp2=(Fp2<G)*1.0
        Fp=Fp1+Fp2 
        Fp=(Fp>1)*1.0 # 取极大值，并置零其他值 
        G_NMS=G*Fp # 保留极大值
        G_NMS=G_NMS.numpy()
        return G_NMS

    # --------- 构建图割 -----------
    def Graphedge(self,img_data):
        '''
        将像素与周围最临近像素合并
        img_data 为梯度值，
        梯度方向是否加入其中
        '''
        # step 1 创建一个图结构
        graph=GraphBase.GenerateNetwork(img_data) # 生成对应的图结构
        # step 2 根据邻接性，度量与周围点的临近性
        graph_near_arr=np.ones((graph.node_num,4),dtype=np.float)*-1 # 确定最小树，这里保留最邻近性，这里保留最临近的点集
        # 开始检索
        # 使用矩阵写法，提升速度
        # 定义九邻域
        img_data_ex=np.pad(img_data,pad_width=1,mode='edge') # 边界外扩展
        # 构建8邻域的边界点
        node_data=np.concatenate([img_data_ex[:-2,:-2].reshape(-1,1),
                                    img_data_ex[:-2,1:-1].reshape(-1,1),
                                    img_data_ex[:-2,2:].reshape(-1,1),
                                    img_data_ex[1:-1,:-2].reshape(-1,1),
                                    img_data_ex[1:-1,1:-1].reshape(-1,1),
                                    img_data_ex[1:-1,2:].reshape(-1,1),
                                    img_data_ex[2:,:-2].reshape(-1,1) ,
                                    img_data_ex[2:,1:-1].reshape(-1,1),
                                    img_data_ex[2:,2:].reshape(-1,1)],axis=1)
        img_ids_ex=np.pad(graph.node_id_arr,pad_width=1,mode='edge') # 边界外扩展
        node_ids=np.concatenate([img_ids_ex[:-2,:-2].reshape(-1,1),
                                    img_ids_ex[:-2,1:-1].reshape(-1,1),
                                    img_ids_ex[:-2,2:].reshape(-1,1),
                                    img_ids_ex[1:-1,:-2].reshape(-1,1),
                                    img_ids_ex[1:-1,1:-1].reshape(-1,1),
                                    img_ids_ex[1:-1,2:].reshape(-1,1),
                                    img_ids_ex[2:,:-2].reshape(-1,1) ,
                                    img_ids_ex[2:,1:-1].reshape(-1,1),
                                    img_ids_ex[2:,2:].reshape(-1,1)],axis=1)
        # 定义临近掩膜版
        max_value=np.max(node_data)-np.min(node_data)+2
        node_mask=np.copy(node_ids)
        node_mask=node_mask==graph.node_id_arr.reshape(-1,1) # 判断与中间
        node_mask=node_mask*max_value
        node_data=node_data+node_mask
        # 定义计算结果
        node_data=node_data-img_data.reshape(-1,1)
        node_data=np.abs(node_data) # 定义运行差值
        node_data_min=np.min(node_data,axis=1).reshape(-1,1)
        node_data=node_data==node_data_min
        node_ids_min=node_ids*node_data
        node_ids_min=np.max(node_ids_min,axis=1).reshape(-1,1) # 选择最有可能的id
        # 保留定义临近值
        graph_near_arr[:,0]=graph.node_id_arr.reshape(-1) # 起点
        graph_near_arr[:,1]=node_ids_min.reshape(-1) # 终点
        graph_near_arr[:,2]=node_data_min.reshape(-1) # 绝对值
        # 最小生成树
        sub_graph_id=0# 定义子类
        # 开始构建连接，并将已经连接的节点记录修改为特定的子类
        # 定义end 结构，方便快速查询
        end_start_dict={}
        for i in range(graph_near_arr.shape[0]):
            start_node=int(graph_near_arr[i,0])
            end_node=int(graph_near_arr[i,1])
            if end_node not in end_start_dict:
                end_start_dict[end_node]=[]
            end_start_dict[end_node].append(start_node)

        for i in range(graph_near_arr.shape[0]):
            if graph_near_arr[i,3]==-1:
                #graph_near_arr=sub_graph_id
                #sub_graph_id=sub_graph_id+1
                # 扩散传播 这里使用队列的方式
                # 向前传播 向后传播
                forward_ids=[graph_near_arr[i,0]]
                next_ids=[graph_near_arr[i,1]] # 向后传播
                while len(forward_ids)>0 or len(next_ids)>0: # 表示当前节点已经处理过了
                    forward_id=int(forward_ids.pop(0)) if len(forward_ids)>0 else None
                    next_id=int(next_ids.pop(0)) if len(next_ids)>0 else None
                    if not  forward_id  is None and graph_near_arr[forward_id,3]!=sub_graph_id: # 向前传播
                        graph_near_arr[forward_id,3]=sub_graph_id
                        forward_id_indexs=end_start_dict[forward_id] if forward_id in end_start_dict else []# 查询速度
                        forward_ids=forward_ids+graph_near_arr[forward_id_indexs,0].reshape(-1).tolist() 
                        next_ids.append(graph_near_arr[forward_id,1]) # 向后传播
                    if not next_id  is None and graph_near_arr[next_id,3]!=sub_graph_id: # 向后传播
                        graph_near_arr[next_id,3]=sub_graph_id
                        next_ids.append(graph_near_arr[next_id,1])   
                        # 定义再次前传
                        forward_id_indexs=end_start_dict[next_id] if next_id in end_start_dict else []
                        forward_ids=forward_ids+graph_near_arr[forward_id_indexs,0].reshape(-1).tolist() 

                sub_graph_id=sub_graph_id+1
            else:
                continue
        # 输出图像
        result=np.copy(graph.node_id_arr)
        result=graph_near_arr[:,3].reshape(graph.node_id_arr.shape[0],graph.node_id_arr.shape[1])
        # 定义结果
        sub_graph_ids=graph_near_arr[:,3].reshape(-1).tolist()
        sub_graph_ids=list(set(sub_graph_ids))
        result_mat=np.copy(result)
        # 构建映射
        ids_ref={}
        for i in sub_graph_ids.copy():
            ids_ref[i]=sub_graph_ids[random.randint(0,len(sub_graph_ids)-1)]

        for i in range(result_mat.shape[0]):
            for j in range(result_mat.shape[1]):
                result_mat[i,j]=ids_ref[result_mat[i,j]]
        scio.savemat('/media/gis/databackup/ayc/modellist/kernel_result/connect.mat',{'connect':result_mat})
        return result 

    def OTSU(self,G_NMS):
        '''
        OTSU 方法自适应滤波,
        做个简单的论证，
        考虑到的场景所使用的影像图像中的，边缘的像素<<内容像素量。
        因此整体来说的这个是属于gamma分布。
        gamma分布可以推导出，其k值分割点，一定在gamma峰值右边，因此牛顿迭代法起点在3/4
        fit_alpha, fit_loc, fit_beta=stats.gamma.fit(data)
        '''
        min_value=np.min(G_NMS)
        max_value=np.max(G_NMS)
        bins=1001
        hist,bin_edges=np.histogram(G_NMS,bins=int(bins))
        bins=list(range(len(bin_edges)))
        bin_edges_dd=np.array(bins[1:])-np.array(bins[0:-1])
        bin_edges_dd=bin_edges_dd/2
        bin_mid_value=bins[0:-1]+bin_edges_dd # 中心点
        

        # 计算mg

        hist_predict=hist/np.sum(hist) # 对应概率
        m=hist_predict*bin_mid_value
        mg=np.sum(hist_predict*bin_mid_value)
        xn=0
        y=-100
        for i in range(len(hist)):
            pk_n=np.sum(hist_predict[:i])
            mk_n=np.sum(m[:i])
            yn=((pk_n*mg-mk_n)**2)/(pk_n*(1-pk_n))
            if y<yn:
                y=yn
                xn=i
        k=bin_edges[xn]
        return k
  
    # 算法主流程
    def computer(self,img):
        img_gray=self.pred_process(img) if not self.pred_process  is None else img[:,:,0]
        h,w=img_gray.shape
        smooth_result=self.smooth_filter(img_gray)
        [G,G_theta]=self.edge_filter(smooth_result)
        #G_NMS=self.NMS(G,G_theta)
        G_NMS=G
        scio.savemat('/media/gis/databackup/ayc/modellist/kernel_result/connect_edeg.mat',{'edge':G_NMS})
        #edge_graph=self.Graphedge(G_NMS)
        k=self.OTSU(G_NMS)
        result=(G_NMS>k)#*G_NMS
        plt.imshow(result)
        plt.show()
        return result

class SAEE(object):
    '''
    基于超像素的结合的算法整理
    '''
    def __init__(self):
        super(SAEE,self).__init__()
        pass

    def computer(img):
        pass

# ------测试模块-----
def Uint_test_canny(img,label,seg):
    '''
    针对图像的中的每个结果进行处理
    这里假设影像数据具有3个特征通道
    '''
    h,w,c=img.shape
    fig=plt.figure(figsize=(80,100))
    img_list=[]
    smooth_gray_ls=[]
    G_ls=[]
    G_NMS_ls=[]
    for c_ix in range(c):
        # 计算流程
        img_gray=img[:,:,c_ix]
        smooth_gray=canny.wmcfilter(img_gray)
        [G,G_theta]=canny.sobelfilter(smooth_gray)
        G_NMS=canny.NMS(G,G_theta)

        smooth_gray_ls.append(smooth_gray)
        G_ls.append(G)
        G_NMS_ls.append(G_NMS)

        plt.subplot(c,6,1+c_ix*6)
        plt.imshow(img)
        plt.title('ori_img')
        plt.axis('off')

        plt.subplot(c,6,2+c_ix*6)
        plt.imshow(img_gray)
        plt.title('img_gray')
        plt.axis('off')

        plt.subplot(c,6,3+c_ix*6)
        plt.imshow(smooth_gray)
        plt.title('smooth_gray')
        plt.axis('off')

        plt.subplot(c,6,4+c_ix*6)
        plt.imshow(G)
        plt.title('G')
        plt.axis('off')

        plt.subplot(c,6,5+c_ix*6)
        plt.imshow(G_NMS)
        plt.title('G_nms')
        plt.axis('off')
        
        plt.subplot(c,6,6+c_ix*6)
        plt.imshow(G_theta)
        plt.title('G_theta')
        plt.axis('off')    
    plt.tight_layout()
    fig=plt.figure(figsize=(80,100)) # 展示每张图的label
    plt.subplot(2,3,1)
    plt.imshow(img)
    plt.axis('off')
    plt.subplot(2,3,2)
    plt.imshow(label)
    plt.axis('off')
    plt.subplot(2,3,3)
    plt.imshow(seg)
    plt.axis('off')
    # 构建    # 边界滤波情况,选取3类型
    plt.subplot(2,3,4)
    edge_1=np.transpose( np.array([label*100,G_NMS_ls[0],G_NMS_ls[0]]),(1,2,0)).astype(np.uint8)
    plt.imshow(edge_1)
    plt.axis('off')   
    plt.subplot(2,3,5)
    edge_2=np.transpose( np.array([label*100,G_NMS_ls[1],G_NMS_ls[1]]),(1,2,0)).astype(np.uint8)
    plt.imshow(edge_1)
    plt.axis('off')   
    plt.subplot(2,3,6)
    edge_2=np.transpose( np.array([label*100,G_NMS_ls[1],G_NMS_ls[1]]),(1,2,0)).astype(np.uint8)
    plt.imshow(edge_1)
    plt.axis('off') 
    plt.tight_layout()
    plt.show()
    # 处理流程



if __name__=="__main__":
    # 单元测试 1445 2037  874
    img=np.load('/media/gis/databackup/ayc/modellist/dataset/nanchang/E512/ALLBands/train/img/1445.npy')
    label=np.load('/media/gis/databackup/ayc/modellist/dataset/nanchang/E512/ALLBands/train/label/1445.npy')
    seg=np.load('/media/gis/databackup/ayc/modellist/dataset/nanchang/E512/ALLBands/train/seg/1445.npy')
    #img=img[:,:,:]
    #Uint_test_canny(img,label,seg)
    cannyFilter=canny(None,smooth_flter="wmc",edge_filter="sobel")
    img_canny=cannyFilter.computer(img)
    pass